# Citizen forecasts of Mexican presidential elections, 2000-2024
# Andreas Murr
# Creates Tables 4-9
# Writes functions invz(), cor.pos(), cor.pos.mul()
# Calls hlm.r

# clear working memory

rm(list = ls())

# write functions

## summary statistics
describe = function(x){
  out = t(apply(x, 2, function(y){c(mean(y), sd(y), min(y), max(y))}))
  colnames(out) = c("Mean", "SD", "Min", "Max") 
  noquote(format(round(out, 2), digits = 2))
}

## correlations

# the inverse of the z transform
invz = function(z){
  # Lindley, D. V., and Scott, W. F., New Cambridge Elementary Statistical Tables, Cambridge: Cambridge University Press (1995) [1st edn (1984)].  Table 17
  # Neave, H. R., Statistics Tables for Mathematicians, Engineers, Economists and the Behavioural and Management Sciences, London: George Allen & Unwin (1978).  Table 6.3
  r = (exp(2 * z) - 1) / (exp(2 * z) + 1)
  return(r)
}

# the posterior of bivariate correlation
cor.pos = function(x, y, alpha = 0.05){
  n = length(x)
  r = cor(x, y)
  z = atanh(r)
  z + c(-1, 1) * qnorm(1 - alpha / 2) * sqrt(1 / n)
  invz(z + c(-1, 1) * qnorm(1 - alpha / 2) * sqrt(1 / n))
}

# posterior of multivariate correlations
cor.pos.mul = function(X, alpha = .05, interval = TRUE){
  # compute correlation
  R = cor(X)
  # create matrix for output
  k = ncol(X)
  O = matrix("", nrow = (k - 1) * 2, ncol = k - 1)
  # compute correlation and interval (standard error)
  i = 1
  j = i + 1
  while (i < k){
    while (j <= k){
      X.new = na.omit(X[,c(i,j)])
      x = X.new[,1]
      y = X.new[,2]
      r = cor(x, y)
      int = cor.pos(x, y, alpha = alpha)
      print.r = formatC(r, format = "f", digits = 2)
      if (interval == FALSE){
        se = (int[2] - int[1]) / (qnorm(1 - alpha / 2) * 2)
        print.i = paste("(", formatC(se, format = "f", digits = 2), ")", sep = "")        
      } else{
        print.i = paste("[", paste(formatC(int, format = "f", digits = 2), collapse = "; "), "]", sep = "")
      }
      # lower triangle
      O[(j-1) + (j-1) - 1, i] = print.r 
      O[(j-1) + (j-1),i] = print.i
      j = j + 1
    }
    i = i + 1
    j = i + 1
  }
  # return results
  noquote(O)
}

# =============
# = load data =
# =============

load("study-3-data.rdata")

# ===========
# = table 4 =
# ===========

d00 = describe(sel00[,-c(1:3,13)])
d06 = describe(sel06[,-c(1:3,13)])
d12 = describe(sel12[,-c(1:3,13)])
d18 = describe(sel18[,-c(1:3,13)])

D = matrix(NA, nrow = nrow(d00)*4, ncol = 6)
for (i in 1:nrow(d00)){
  D[(1+(i-1)*4),1] = rownames(d00)[i]
  D[(1+(i-1)*4):(i*4),2] = c("2000", "2006", "2012", "2018")
  D[(1+(i-1)*4):(i*4),3:6] = rbind(d00[i,], d06[i,], d12[i,], d18[i,])
}
noquote(D)

# ===========
# = table 5 =
# ===========

C00 = sel00[,-c(1:3,5)]
C06 = sel06[,-c(1:3,5)]
C12 = sel12[,-c(1:3,5)]
C00$pid.mrn = NA
C06$pid.mrn = NA
C12$pid.mrn = NA
C18 = cbind(sel18[,c(4,6:8)], NA, sel18[,c(10:13, 9)])
names(C18)[5] = "pid.prd"
C = rbind(C00, C06, C12, C18)
C = C[,c(1,9,2:3,10,4:6,7:8)]

C.out = cor.pos.mul(C, interval = FALSE)
colnames(C.out) = colnames(C)[-ncol(C)]
rownames(C.out) = rep("", nrow(C.out))
rownames(C.out)[seq(1, nrow(C.out), 2)] = colnames(C)[-1]
C.out

# =============================================
# = fit bayesian heteroskedastic linear model =
# =============================================

# create model matrix

y00.pan = sel00$err.pan
y06.pan = sel06$err.pan
y12.pan = sel12$err.pan
y18.pan = sel18$err.pan

y00.prd = sel00$err.prd
y06.prd = sel06$err.prd
y12.prd = sel12$err.prd
y18.mrn = sel18$err.mrn

y00.pri = sel00$err.pri
y06.pri = sel06$err.pri
y12.pri = sel12$err.pri
y18.pri = sel18$err.pri

X00.pan = with(sel00, cbind(1, female, age100, educ, interest, pid.pan, clean2, time))
X06.pan = with(sel06, cbind(1, female, age100, educ, interest, pid.pan, clean2, time))
X12.pan = with(sel12, cbind(1, female, age100, educ, interest, pid.pan, clean2, time))
X18.pan = with(sel18, cbind(1, female, age100, educ, interest, pid.pan, clean2, time))

X00.prd = with(sel00, cbind(1, female, age100, educ, interest, pid.prd, clean2, time))
X06.prd = with(sel06, cbind(1, female, age100, educ, interest, pid.prd, clean2, time))
X12.prd = with(sel12, cbind(1, female, age100, educ, interest, pid.prd, clean2, time))
X18.mrn = with(sel18, cbind(1, female, age100, educ, interest, pid.mrn, clean2, time))

X00.pri = with(sel00, cbind(1, female, age100, educ, interest, pid.pri, clean2, time))
X06.pri = with(sel06, cbind(1, female, age100, educ, interest, pid.pri, clean2, time))
X12.pri = with(sel12, cbind(1, female, age100, educ, interest, pid.pri, clean2, time))
X18.pri = with(sel18, cbind(1, female, age100, educ, interest, pid.pri, clean2, time))

Z00.pan = X00.pan
Z06.pan = X06.pan
Z12.pan = X12.pan
Z18.pan = X18.pan

Z00.prd = X00.prd
Z06.prd = X06.prd
Z12.prd = X12.prd
Z18.mrn = X18.mrn

Z00.pri = X00.pri
Z06.pri = X06.pri
Z12.pri = X12.pri
Z18.pri = X18.pri

# compute heteroskedastic linear model

source("hlm.r")

## fit models

  set.seed(512438)
  # pan
  fit00.pan = hlm(y00.pan, X00.pan, Z00.pan, n.iter = 1000, n.prog = 100)
  fit06.pan = hlm(y06.pan, X06.pan, Z06.pan, n.iter = 1000, n.prog = 100)
  fit12.pan = hlm(y12.pan, X12.pan, Z12.pan, n.iter = 1000, n.prog = 100)
  fit18.pan = hlm(y18.pan, X18.pan, Z18.pan, n.iter = 1000, n.prog = 100)
  # prd / mrn
  fit00.prd = hlm(y00.prd, X00.prd, Z00.prd, n.iter = 1000, n.prog = 100)
  fit06.prd = hlm(y06.prd, X06.prd, Z06.prd, n.iter = 1000, n.prog = 100)
  fit12.prd = hlm(y12.prd, X12.prd, Z12.prd, n.iter = 1000, n.prog = 100)
  fit18.mrn = hlm(y18.mrn, X18.mrn, Z18.mrn, n.iter = 1000, n.prog = 100)
  # pri
  fit00.pri = hlm(y00.pri, X00.pri, Z00.pri, n.iter = 1000, n.prog = 100)
  fit06.pri = hlm(y06.pri, X06.pri, Z06.pri, n.iter = 1000, n.prog = 100)
  fit12.pri = hlm(y12.pri, X12.pri, Z12.pri, n.iter = 1000, n.prog = 100)
  fit18.pri = hlm(y18.pri, X18.pri, Z18.pri, n.iter = 1000, n.prog = 100)

## check convergence

  # pan
  plot.hlm.iter(fit00.pan)
  plot.hlm.iter(fit00.pan, n.burn = 100, store = 1)
  plot.hlm.iter(fit06.pan)
  plot.hlm.iter(fit06.pan, n.burn = 100, store = 1)
  plot.hlm.iter(fit12.pan)
  plot.hlm.iter(fit12.pan, n.burn = 100, store = 1)
  plot.hlm.iter(fit18.pan)
  plot.hlm.iter(fit18.pan, n.burn = 100, store = 1)
  # prd / mrn
  plot.hlm.iter(fit00.prd)
  plot.hlm.iter(fit00.prd, n.burn = 100, store = 1)
  plot.hlm.iter(fit06.prd)
  plot.hlm.iter(fit06.prd, n.burn = 100, store = 1)
  plot.hlm.iter(fit12.prd)
  plot.hlm.iter(fit12.prd, n.burn = 100, store = 1)
  plot.hlm.iter(fit18.mrn)
  plot.hlm.iter(fit18.mrn, n.burn = 100, store = 1)
  # pri
  plot.hlm.iter(fit00.pri)
  plot.hlm.iter(fit00.pri, n.burn = 100, store = 1)
  plot.hlm.iter(fit06.pri)
  plot.hlm.iter(fit06.pri, n.burn = 100, store = 1)
  plot.hlm.iter(fit12.pri)
  plot.hlm.iter(fit12.pri, n.burn = 100, store = 1)
  plot.hlm.iter(fit18.pri)
  plot.hlm.iter(fit18.pri, n.burn = 100, store = 1)

# ===========
# = table 6 =
# ===========

r = matrix(NA, nrow = (ncol(X00.pan) - 1) * 2, ncol = 3)
rownames(r) = rep("", nrow(r))
j = 0 
for (i in 2:ncol(X00.pan)){
  # create range the variable
  r00 = range(X00.pan[,i])
  r06 = range(X06.pan[,i])
  r12 = range(X12.pan[,i])
  r18 = range(X18.pan[,i])
  lo = max(c(r00[1], r06[1], r12[1], r18[1]))
  hi = min(c(r00[2], r06[2], r12[2], r18[2]))
  # create two copies of data
  ## pan
  X00.pan.lo = X00.pan.hi = X00.pan
  X06.pan.lo = X06.pan.hi = X06.pan
  X12.pan.lo = X12.pan.hi = X12.pan
  X18.pan.lo = X18.pan.hi = X18.pan
  ## prd / mrn
  X00.prd.lo = X00.prd.hi = X00.prd
  X06.prd.lo = X06.prd.hi = X06.prd
  X12.prd.lo = X12.prd.hi = X12.prd
  X18.mrn.lo = X18.mrn.hi = X18.mrn
  ## pri
  X00.pri.lo = X00.pri.hi = X00.pri
  X06.pri.lo = X06.pri.hi = X06.pri
  X12.pri.lo = X12.pri.hi = X12.pri
  X18.pri.lo = X18.pri.hi = X18.pri
  # set predictor to low value
  X00.pan.lo[,i] = lo
  X06.pan.lo[,i] = lo
  X12.pan.lo[,i] = lo
  X18.pan.lo[,i] = lo
  X00.prd.lo[,i] = lo
  X06.prd.lo[,i] = lo
  X12.prd.lo[,i] = lo
  X18.mrn.lo[,i] = lo
  X00.pri.lo[,i] = lo
  X06.pri.lo[,i] = lo
  X12.pri.lo[,i] = lo
  X18.pri.lo[,i] = lo
  # set predictor to high value 
  X00.pan.hi[,i] = hi
  X06.pan.hi[,i] = hi
  X12.pan.hi[,i] = hi
  X18.pan.hi[,i] = hi
  X00.prd.hi[,i] = hi
  X06.prd.hi[,i] = hi
  X12.prd.hi[,i] = hi
  X18.mrn.hi[,i] = hi
  X00.pri.hi[,i] = hi
  X06.pri.hi[,i] = hi
  X12.pri.hi[,i] = hi
  X18.pri.hi[,i] = hi
  # compute bias squared
  ## for low values
  bias.sq.00.pan.lo = (X00.pan.lo%*%t(fit00.pan[[1]]))^2
  bias.sq.06.pan.lo = (X06.pan.lo%*%t(fit06.pan[[1]]))^2
  bias.sq.12.pan.lo = (X12.pan.lo%*%t(fit12.pan[[1]]))^2
  bias.sq.18.pan.lo = (X18.pan.lo%*%t(fit18.pan[[1]]))^2
  bias.sq.00.prd.lo = (X00.prd.lo%*%t(fit00.prd[[1]]))^2
  bias.sq.06.prd.lo = (X06.prd.lo%*%t(fit06.prd[[1]]))^2
  bias.sq.12.prd.lo = (X12.prd.lo%*%t(fit12.prd[[1]]))^2
  bias.sq.18.mrn.lo = (X18.mrn.lo%*%t(fit18.mrn[[1]]))^2
  bias.sq.00.pri.lo = (X00.pri.lo%*%t(fit00.pri[[1]]))^2
  bias.sq.06.pri.lo = (X06.pri.lo%*%t(fit06.pri[[1]]))^2
  bias.sq.12.pri.lo = (X12.pri.lo%*%t(fit12.pri[[1]]))^2
  bias.sq.18.pri.lo = (X18.pri.lo%*%t(fit18.pri[[1]]))^2
  ## for high values
  bias.sq.00.pan.hi = (X00.pan.hi%*%t(fit00.pan[[1]]))^2
  bias.sq.06.pan.hi = (X06.pan.hi%*%t(fit06.pan[[1]]))^2
  bias.sq.12.pan.hi = (X12.pan.hi%*%t(fit12.pan[[1]]))^2
  bias.sq.18.pan.hi = (X18.pan.hi%*%t(fit18.pan[[1]]))^2
  bias.sq.00.prd.hi = (X00.prd.hi%*%t(fit00.prd[[1]]))^2
  bias.sq.06.prd.hi = (X06.prd.hi%*%t(fit06.prd[[1]]))^2
  bias.sq.12.prd.hi = (X12.prd.hi%*%t(fit12.prd[[1]]))^2
  bias.sq.18.mrn.hi = (X18.mrn.hi%*%t(fit18.mrn[[1]]))^2
  bias.sq.00.pri.hi = (X00.pri.hi%*%t(fit00.pri[[1]]))^2
  bias.sq.06.pri.hi = (X06.pri.hi%*%t(fit06.pri[[1]]))^2
  bias.sq.12.pri.hi = (X12.pri.hi%*%t(fit12.pri[[1]]))^2
  bias.sq.18.pri.hi = (X18.pri.hi%*%t(fit18.pri[[1]]))^2
  ## compute difference
  d.bias.sq.00.pan = bias.sq.00.pan.hi - bias.sq.00.pan.lo
  d.bias.sq.06.pan = bias.sq.06.pan.hi - bias.sq.06.pan.lo
  d.bias.sq.12.pan = bias.sq.12.pan.hi - bias.sq.12.pan.lo
  d.bias.sq.18.pan = bias.sq.18.pan.hi - bias.sq.18.pan.lo
  d.bias.sq.00.prd = bias.sq.00.prd.hi - bias.sq.00.prd.lo
  d.bias.sq.06.prd = bias.sq.06.prd.hi - bias.sq.06.prd.lo
  d.bias.sq.12.prd = bias.sq.12.prd.hi - bias.sq.12.prd.lo
  d.bias.sq.18.mrn = bias.sq.18.mrn.hi - bias.sq.18.mrn.lo
  d.bias.sq.00.pri = bias.sq.00.pri.hi - bias.sq.00.pri.lo
  d.bias.sq.06.pri = bias.sq.06.pri.hi - bias.sq.06.pri.lo
  d.bias.sq.12.pri = bias.sq.12.pri.hi - bias.sq.12.pri.lo
  d.bias.sq.18.pri = bias.sq.18.pri.hi - bias.sq.18.pri.lo
  ## compute individual mean in each election for each parties
  m01 = mean(apply(d.bias.sq.00.pan, 1, mean))
  m02 = mean(apply(d.bias.sq.06.pan, 1, mean))
  m03 = mean(apply(d.bias.sq.12.pan, 1, mean))
  m04 = mean(apply(d.bias.sq.18.pan, 1, mean))
  m05 = mean(apply(d.bias.sq.00.prd, 1, mean))
  m06 = mean(apply(d.bias.sq.06.prd, 1, mean))
  m07 = mean(apply(d.bias.sq.12.prd, 1, mean))
  m08 = mean(apply(d.bias.sq.18.mrn, 1, mean))
  m09 = mean(apply(d.bias.sq.00.prd, 1, mean))
  m10 = mean(apply(d.bias.sq.06.pri, 1, mean))
  m11 = mean(apply(d.bias.sq.12.pri, 1, mean))
  m12 = mean(apply(d.bias.sq.18.pri, 1, mean))
  ## compute variance of mean 
  ## var(x.bar) = var(1/n sum_i^n x_i) 
  ## = 1/n^2 sum_i^n var(x_i) = 1 / n^2 sum_i^n sigma_i^2
  v01 = sum(apply(d.bias.sq.00.pan, 1, var)) / nrow(d.bias.sq.00.pan)^2
  v02 = sum(apply(d.bias.sq.06.pan, 1, var)) / nrow(d.bias.sq.06.pan)^2
  v03 = sum(apply(d.bias.sq.12.pan, 1, var)) / nrow(d.bias.sq.12.pan)^2
  v04 = sum(apply(d.bias.sq.18.pan, 1, var)) / nrow(d.bias.sq.18.pan)^2
  v05 = sum(apply(d.bias.sq.00.prd, 1, var)) / nrow(d.bias.sq.00.prd)^2
  v06 = sum(apply(d.bias.sq.06.prd, 1, var)) / nrow(d.bias.sq.06.prd)^2
  v07 = sum(apply(d.bias.sq.12.prd, 1, var)) / nrow(d.bias.sq.12.prd)^2
  v08 = sum(apply(d.bias.sq.18.mrn, 1, var)) / nrow(d.bias.sq.18.mrn)^2
  v09 = sum(apply(d.bias.sq.00.prd, 1, var)) / nrow(d.bias.sq.00.prd)^2
  v10 = sum(apply(d.bias.sq.06.pri, 1, var)) / nrow(d.bias.sq.06.pri)^2
  v11 = sum(apply(d.bias.sq.12.pri, 1, var)) / nrow(d.bias.sq.12.pri)^2
  v12 = sum(apply(d.bias.sq.18.pri, 1, var)) / nrow(d.bias.sq.18.pri)^2
  ## save results
  bias.sq.m = (m01 + m02 + m03 + m04 + m05 + m06 + m07 + m08 + m09 + m10 + m11 + m12) / 12
  bias.sq.s = sqrt((v01 + v02 + v03 + v04 + v05 + v06 + v07 + v08 + v09 + v10 + v11 + v12) / 12^2)
  # compute variance
  ## for low values
  var.00.pan.lo = exp(X00.pan.lo%*%t(fit00.pan[[2]]))
  var.06.pan.lo = exp(X06.pan.lo%*%t(fit06.pan[[2]]))
  var.12.pan.lo = exp(X12.pan.lo%*%t(fit12.pan[[2]]))
  var.18.pan.lo = exp(X18.pan.lo%*%t(fit18.pan[[2]]))
  var.00.prd.lo = exp(X00.prd.lo%*%t(fit00.prd[[2]]))
  var.06.prd.lo = exp(X06.prd.lo%*%t(fit06.prd[[2]]))
  var.12.prd.lo = exp(X12.prd.lo%*%t(fit12.prd[[2]]))
  var.18.mrn.lo = exp(X18.mrn.lo%*%t(fit18.mrn[[2]]))
  var.00.pri.lo = exp(X00.pri.lo%*%t(fit00.pri[[2]]))
  var.06.pri.lo = exp(X06.pri.lo%*%t(fit06.pri[[2]]))
  var.12.pri.lo = exp(X12.pri.lo%*%t(fit12.pri[[2]]))
  var.18.pri.lo = exp(X18.pri.lo%*%t(fit18.pri[[2]]))
  ## for high values
  var.00.pan.hi = exp(X00.pan.hi%*%t(fit00.pan[[2]]))
  var.06.pan.hi = exp(X06.pan.hi%*%t(fit06.pan[[2]]))
  var.12.pan.hi = exp(X12.pan.hi%*%t(fit12.pan[[2]]))
  var.18.pan.hi = exp(X18.pan.hi%*%t(fit18.pan[[2]]))
  var.00.prd.hi = exp(X00.prd.hi%*%t(fit00.prd[[2]]))
  var.06.prd.hi = exp(X06.prd.hi%*%t(fit06.prd[[2]]))
  var.12.prd.hi = exp(X12.prd.hi%*%t(fit12.prd[[2]]))
  var.18.mrn.hi = exp(X18.mrn.hi%*%t(fit18.mrn[[2]]))
  var.00.pri.hi = exp(X00.pri.hi%*%t(fit00.pri[[2]]))
  var.06.pri.hi = exp(X06.pri.hi%*%t(fit06.pri[[2]]))
  var.12.pri.hi = exp(X12.pri.hi%*%t(fit12.pri[[2]]))
  var.18.pri.hi = exp(X18.pri.hi%*%t(fit18.pri[[2]]))
  ## compute difference
  d.var.00.pan = var.00.pan.hi - var.00.pan.lo
  d.var.06.pan = var.06.pan.hi - var.06.pan.lo
  d.var.12.pan = var.12.pan.hi - var.12.pan.lo
  d.var.18.pan = var.18.pan.hi - var.18.pan.lo
  d.var.00.prd = var.00.prd.hi - var.00.prd.lo
  d.var.06.prd = var.06.prd.hi - var.06.prd.lo
  d.var.12.prd = var.12.prd.hi - var.12.prd.lo
  d.var.18.mrn = var.18.mrn.hi - var.18.mrn.lo
  d.var.00.pri = var.00.pri.hi - var.00.pri.lo
  d.var.06.pri = var.06.pri.hi - var.06.pri.lo
  d.var.12.pri = var.12.pri.hi - var.12.pri.lo
  d.var.18.pri = var.18.pri.hi - var.18.pri.lo
  ## compute individual mean in each election for each parties
  m01 = mean(apply(d.var.00.pan, 1, mean))
  m02 = mean(apply(d.var.06.pan, 1, mean))
  m03 = mean(apply(d.var.12.pan, 1, mean))
  m04 = mean(apply(d.var.18.pan, 1, mean))
  m05 = mean(apply(d.var.00.prd, 1, mean))
  m06 = mean(apply(d.var.06.prd, 1, mean))
  m07 = mean(apply(d.var.12.prd, 1, mean))
  m08 = mean(apply(d.var.18.mrn, 1, mean))
  m09 = mean(apply(d.var.00.prd, 1, mean))
  m10 = mean(apply(d.var.06.pri, 1, mean))
  m11 = mean(apply(d.var.12.pri, 1, mean))
  m12 = mean(apply(d.var.18.pri, 1, mean))
  ## compute variance of mean 
  ## var(x.bar) = var(1/n sum_i^n x_i) 
  ## = 1/n^2 sum_i^n var(x_i) = 1 / n^2 sum_i^n sigma_i^2
  v01 = sum(apply(d.var.00.pan, 1, var)) / nrow(d.var.00.pan)^2
  v02 = sum(apply(d.var.06.pan, 1, var)) / nrow(d.var.06.pan)^2
  v03 = sum(apply(d.var.12.pan, 1, var)) / nrow(d.var.12.pan)^2
  v04 = sum(apply(d.var.18.pan, 1, var)) / nrow(d.var.18.pan)^2
  v05 = sum(apply(d.var.00.prd, 1, var)) / nrow(d.var.00.prd)^2
  v06 = sum(apply(d.var.06.prd, 1, var)) / nrow(d.var.06.prd)^2
  v07 = sum(apply(d.var.12.prd, 1, var)) / nrow(d.var.12.prd)^2
  v08 = sum(apply(d.var.18.mrn, 1, var)) / nrow(d.var.18.mrn)^2
  v09 = sum(apply(d.var.00.prd, 1, var)) / nrow(d.var.00.prd)^2
  v10 = sum(apply(d.var.06.pri, 1, var)) / nrow(d.var.06.pri)^2
  v11 = sum(apply(d.var.12.pri, 1, var)) / nrow(d.var.12.pri)^2
  v12 = sum(apply(d.var.18.pri, 1, var)) / nrow(d.var.18.pri)^2
  ## save results
  var.sq.m = (m01 + m02 + m03 + m04 + m05 + m06 + m07 + m08 + m09 + m10 + m11 + m12) / 12
  var.sq.s = sqrt((v01 + v02 + v03 + v04 + v05 + v06 + v07 + v08 + v09 + v10 + v11 + v12) / 12^2)
  # compute difference in mean squared error
  ## is difference in bias.sq + difference in variance
  m01 = mean(apply(d.bias.sq.00.pan + d.var.00.pan, 1, mean))
  m02 = mean(apply(d.bias.sq.06.pan + d.var.06.pan, 1, mean))
  m03 = mean(apply(d.bias.sq.12.pan + d.var.12.pan, 1, mean))
  m04 = mean(apply(d.bias.sq.18.pan + d.var.18.pan, 1, mean))
  m05 = mean(apply(d.bias.sq.00.prd + d.var.00.prd, 1, mean))
  m06 = mean(apply(d.bias.sq.06.prd + d.var.06.prd, 1, mean))
  m07 = mean(apply(d.bias.sq.12.prd + d.var.12.prd, 1, mean))
  m08 = mean(apply(d.bias.sq.18.mrn + d.var.18.mrn, 1, mean))
  m09 = mean(apply(d.bias.sq.00.prd + d.var.00.prd, 1, mean))
  m10 = mean(apply(d.bias.sq.06.pri + d.var.06.pri, 1, mean))
  m11 = mean(apply(d.bias.sq.12.pri + d.var.12.pri, 1, mean))
  m12 = mean(apply(d.bias.sq.18.pri + d.var.18.pri, 1, mean))
  ## compute variance
  v01 = sum(apply(d.bias.sq.00.pan + d.var.00.pan, 1, var)) / nrow(d.bias.sq.00.pan)^2
  v02 = sum(apply(d.bias.sq.06.pan + d.var.06.pan, 1, var)) / nrow(d.bias.sq.06.pan)^2
  v03 = sum(apply(d.bias.sq.12.pan + d.var.12.pan, 1, var)) / nrow(d.bias.sq.12.pan)^2
  v04 = sum(apply(d.bias.sq.18.pan + d.var.18.pan, 1, var)) / nrow(d.bias.sq.18.pan)^2
  v05 = sum(apply(d.bias.sq.00.prd + d.var.00.prd, 1, var)) / nrow(d.bias.sq.00.prd)^2
  v06 = sum(apply(d.bias.sq.06.prd + d.var.06.prd, 1, var)) / nrow(d.bias.sq.06.prd)^2
  v07 = sum(apply(d.bias.sq.12.prd + d.var.12.prd, 1, var)) / nrow(d.bias.sq.12.prd)^2
  v08 = sum(apply(d.bias.sq.18.mrn + d.var.18.mrn, 1, var)) / nrow(d.bias.sq.18.mrn)^2
  v09 = sum(apply(d.bias.sq.00.prd + d.var.00.prd, 1, var)) / nrow(d.bias.sq.00.prd)^2
  v10 = sum(apply(d.bias.sq.06.pri + d.var.06.pri, 1, var)) / nrow(d.bias.sq.06.pri)^2
  v11 = sum(apply(d.bias.sq.12.pri + d.var.12.pri, 1, var)) / nrow(d.bias.sq.12.pri)^2
  v12 = sum(apply(d.bias.sq.18.pri + d.var.18.pri, 1, var)) / nrow(d.bias.sq.18.pri)^2
  ## save results
  mse.sq.m = (m01 + m02 + m03 + m04 + m05 + m06 + m07 + m08 + m09 + m10 + m11 + m12) / 12
  mse.sq.s = sqrt((v01 + v02 + v03 + v04 + v05 + v06 + v07 + v08 + v09 + v10 + v11 + v12) / 12^2)
  ## return results
  r[(i-1) + j,] = c(bias.sq.m, var.sq.m, mse.sq.m)
  r[i + j,] = c(bias.sq.s, var.sq.s, mse.sq.s)
  rownames(r)[(i-1) + j] = paste(colnames(X00.pan)[i], " (", lo, " to ", hi, ")", sep="")
  j = j + 1
}

colnames(r) = c("bias.sq", "variance", "mse")
rownames(r) = rep("", nrow(r))
rownames(r)[seq(1, nrow(r), 2)] = colnames(X00.pan)[-1]

p = r
for (i in seq(1, nrow(p), 2)){
  p[i,] = paste(formatC(r[i,], format = "f", digits = 1))
}
for (i in seq(2, nrow(p), 2)){
  p[i,] = paste("(", formatC(r[i,], format = "f", digits = 1), ")", sep = "")
}
noquote(p)

# return changes in predictor vairables (1st column)

l00 = apply(X00.pan[,-1], 2, min)
h00 = apply(X00.pan[,-1], 2, max)
l06 = apply(X06.pan[,-1], 2, min)
h06 = apply(X06.pan[,-1], 2, max)
l12 = apply(X12.pan[,-1], 2, min)
h12 = apply(X12.pan[,-1], 2, max)
l18 = apply(X18.pan[,-1], 2, min)
h18 = apply(X18.pan[,-1], 2, max)

apply(rbind(l00, l06, l12, l18), 2, max)
apply(rbind(h00, h06, h12, h18), 2, min)

# ==============
# = tables 7-9 =
# ==============

# tabulate summaries of posterior

  # name intercept
  colnames(X00.pan)[1] = "Intercept"
  colnames(X06.pan)[1] = "Intercept"
  colnames(X12.pan)[1] = "Intercept"
  colnames(X18.pan)[1] = "Intercept"
  colnames(X00.prd)[1] = "Intercept"
  colnames(X06.prd)[1] = "Intercept"
  colnames(X12.prd)[1] = "Intercept"
  colnames(X18.mrn)[1] = "Intercept"
  colnames(X00.pri)[1] = "Intercept"
  colnames(X06.pri)[1] = "Intercept"
  colnames(X12.pri)[1] = "Intercept"
  colnames(X18.pri)[1] = "Intercept"

  # table 7
  noquote(cbind(
    summary.hlm(fit00.pan, y00.pan, X00.pan, n.burn = 100, store = 1),
    summary.hlm(fit06.pan, y06.pan, X06.pan, n.burn = 100, store = 1),
    summary.hlm(fit12.pan, y12.pan, X12.pan, n.burn = 100, store = 1),
    summary.hlm(fit18.pan, y18.pan, X18.pan, n.burn = 100, store = 1)
  ))
  # table 8
  noquote(cbind(
    summary.hlm(fit00.prd, y00.prd, X00.prd, n.burn = 100, store = 1),
    summary.hlm(fit06.prd, y06.prd, X06.prd, n.burn = 100, store = 1),
    summary.hlm(fit12.prd, y12.prd, X12.prd, n.burn = 100, store = 1),
    summary.hlm(fit18.mrn, y18.mrn, X18.mrn, n.burn = 100, store = 1)
  ))
  # table 9
  noquote(cbind(
    summary.hlm(fit00.pri, y00.pri, X00.pri, n.burn = 100, store = 1),
    summary.hlm(fit06.pri, y06.pri, X06.pri, n.burn = 100, store = 1),
    summary.hlm(fit12.pri, y12.pri, X12.pri, n.burn = 100, store = 1),
    summary.hlm(fit18.pri, y18.pri, X18.pri, n.burn = 100, store = 1)
  )) 

# ===================
# = end source code =
# ===================